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Abstract 

The development of a gridless computational fluid dynamics (CFD) method for the solution of the 
two-dimensional Euler and Navicr-Stokes equations is described. The method uses only clouds of points 
and does not require that the points be connected to form a grid as is necessary in conventional CFD 
algorithms. The gridless CFD approach appears to resolve the problems and inefficiencies encountered 
with structured or unstructured grid methods, and consequently offers the greatest potential for accurately 
and efficiently solving viscous flows about complex aircraft configurations. The method is described in 
detail and calculations are presented for standard Euler and Navier-Stokes cases to assess the accuracy 
and efficiency of the capability. 


Introduction 

Considerable progress in developing computational fluid dynamics (CFD) methods for aerodynamic 
analysis has been made over the past two decades. 1 The majority of work that has been done in CFD 
over the years has been on developing methods for use on computational grids that have an underlying 
geometrical structure and thus the grids are referred to as “structured”. For example, Fig. 1(a) shows 
a structured grid for the NACA 0012 airfoil. The grid is of C-type topology, has 159 points in the 
wraparound direction, and 49 points in the outward direction. Methods developed for structured grids have 
been applied to a wide variety of geometrical configurations ranging from simple, analytically defined 
airfoil sections such as the NACA 0012 airfoil to complex aircraft configurations such as the F-16A 
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fighter. 2 Although applications of structured grid methods to complex configurations are indeed possible 
they generally require more sophisticated meshing methodologies such as blocked, patched, chimera, 
or hybrid-type grids. For example, the F-16A fighter calculations reported in Ref. 2, which included 
the engine inlet and boundary layer diverter as well as the wing, fuselage, and tail in the geometrical 
modeling, used 27 blocks of structured cells to make up the grid. These more sophisticated meshing 
methodologies, in turn, significantly complicate the solution algorithms of the structured grid methods. 

An alternative approach is the use of unstructured grids. 3-7 In two dimensions, unstructured grids 
typically are constructed from triangles, and in three dimensions, they consist of tetrahedral cells. The 
triangles or tetrahedra may be oriented in an arbitrary way to conform to the geometry, thus making 
it possible to easily generate grids about very complicated shapes. Although not a complicated shape, 
Fig. 1(b) shows an example of an unstructured grid for the NACA 0012 airfoil. The total grid has 
3300 nodes and 6466 triangles. An advantage of methods developed for unstructured grids is that they 
may be applied to complex aircraft configurations without having to make changes to the basic solution 
algorithm. Numerous calculations for complex configurations performed using various Euler codes have 
been reported by several researchers. 3-7 However, applications to three-dimensional configurations using 
unstructured grid Euler codes have tended to be inefficient because the meshes have an excessively large 
number of cells. The excessive number of cells is due, in part, to the current state-of-the-art in generation 
of unstructured tetrahedral grids, which produces meshes that are much finer in the spanwise direction (for 
a given streamwise density) than is necessary for accurate flow computation. To alleviate the problem, 
the cells may be stretched in the spanwise direction when generating the mesh to reduce the number of 
cells. However, the stretching can create convergence and accuracy problems for the flow solver. The 
basic problem is that the tetrahedron is an inefficient geometrical shape (whereas the triangle tends to 
be an efficient shape in two dimensions). A more efficient shape for an isolated wing application is a 
prismatic cell defined by a polyhedron with a triangular cross-section. A mesh of this type uses triangles 
which form prisms when connected in the spanwise direction to grid the planes of the airfoil sections of 
the wing. This approach, though, not only puts structure back into the mesh it is not generally applicable 
to complex three-dimensional configurations. 
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Another problem with the unstructured-grid methodology is encountered in extending the methods for 
solving the Euler equations to the solution of the Navier-Stokes equations, especially in three dimensions. 
For viscous applications, grids generally need to be fine near the body in the outward direction to resolve 
the boundary layer but less fine in the direction along the surface of the body. This naturally leads 
to cells of high aspect ratio which tends to exacerbate the inefficiency of three-dimensional solution 
algorithms based on tetrahcdra. Specifically, the use of tetrahedra for viscous flow applications results 
in an unreasonably large number of cells. The number of cells is in fact absurdly large in comparison to 
grids that are generated for Euler calculations (which are already inefficient because of a large number of 
cells as previously discussed) because of the additional requirement that the mesh be fine near the body. 
To alleviate this problem, a hybrid approach has been developed recently using prismatic cells for the 
solution of the Navier-Stokes equations . 8 In this approach, the surface of the geometry under consideration 
and the outer boundaries of the mesh are gridded using triangles, and instead of generating tetrahedra to 
fill the interior of the computational domain, the triangles on the inner and outer boundaries of the mesh 
arc connected to form prisms. The prisms, of course, require the same number of triangles on the inner 
and outer boundaries. While this hybrid approach is a viable solution to alleviate the inefficiency created 
by using tetrahedral cells to solve the Navier-Stokes equations, it is not necessarily the best approach, 
since it again puts structure back into the mesh and limits some of the advantages of the unstructured 
grid methodology, such as spatial adaption. 

What is truly required to advance the CFD technology to treat complex configurations in viscous 
flows is not to take a step backward toward grid structure, but to take a bold step forward to develop 
methods that do not require the use of grids at all. Hence the solution to the above-mentioned problems 
with structured and unstructured grids is the development of algorithms for solving the Navier-Stokes 
equations based on using only grid points and not on the connectivity information that relates all of the 
points to one another. This type of approach, which may be referred to as “gridless” CFD, has distinct 
advantages over methods that require grids. Since only points are required, or specifically clouds of 
points as suggested by Chakravarthy , 9 gridless CFD methods offer the greatest potential for accurately 
and efficiently solving viscous flows about complex aircraft configurations. It is noted parenthetically, 
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that if finally the grid points too were not required by the solution algorithm, then the ultimate flexibility 
in methodology could be attained. This type of method might then be referred to as “pointless” CFD. 

The purpose of the paper is to report the development of a gridless method for the solution of the 
two-dimensional Euler and Navicr-Stokcs equations. The method uses only clouds of points and does 
not require that the points be connected to form a grid as is necessary in conventional CFD algorithms. 
The governing partial differential equations (PDEs) are solved directly, by performing local least-squares 
curve fits in each cloud of points, and then analytically differentiating the resulting curve-fit equations to 
approximate the derivatives of the PDEs. The method is neither a finite-difference nor a finite-volume 
type approach since differences, metrics, lengths, areas, or volumes are not computed. The method is 
described in further detail and calculations are presented for standard cases to assess the accuracy and 
efficiency of the capability. 


Governing Equations 

In this study the flow is assumed to be governed by the two-dimensional laminar Navier-Stokes 
equations which may be written in differential form as 


£+ k < j -*>+£<'-*>= 0 


where Q is the vector of conserved variables given by 

Q = 

E and F arc the inviscid fluxes in the x and y directions, respectively, defined by 


E= ( 
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F= ( 
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In the viscous fluxes the shear stresses and heat flux terms are defined by 
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In these equations, M^ is the frcestream 
number, and f.i is the molecular viscosity determined using Sutherland’s law. The Euler equations are 
obtained by setting the viscous fluxes equal to zero. 
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Mach number, Re is the Reynolds number, Pr is the Prandtl 



and E v and F v are the viscous fluxes in th 
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Spatial Discretization 


Derivatives 

The spatial derivatives in the governing equations (Eq. (1)) are approximated as follows. In each 
cloud of points, each term of the fluxes is assumed to vary linearly according to 

f(x,y) = a 0 + aix + a 2 y (2) 


where the coefficients ao, a { , and a 2 arc determined from a least-squares curve fit. Performing a least- 
squares fit in a given cloud results in three equations represented in matrix form by 
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where n is the number of points in the cloud and the summations are taken over the n points. The 
solution of Eqs. (3) requires the inversion of a 3 x 3 matrix which is performed for every cloud in the 
computational domain. Having solved these equations for oq , a i , and a 2 , the spatial derivatives are now 
known since by differentiating Eq. (2) it is obvious that 




(4) 


In addition to approximating the spatial derivatives of the governing equations by differentiation of 
the least-squares curve fits, the shear stresses and heat flux terms are calculated the same way. Since 
these terms involve first derivatives of the velocity components or pressure divided by density, the shear 
stresses and heat fluxes can be approximated by defining / to be equal to u, v, or p/p, evaluating the 
terms of Eqs. (3), and inverting the left-hand-side matrix. The resulting values for ai and 02 are the 
derivatives of the specified quantity with respect to x and y, respectively, within a given cloud of points. 

Artificial Dissipation 

The unsteady Euler equations are a set of nondissipative hyperbolic conservation laws that require 
some form of artificial dissipation to prevent oscillations near shock waves and to damp high-frequency 
uncoupled error modes. The unsteady Navicr-Stokes equations also require artificial dissipation for 
similar reasons because the physical viscosity generally is limited to the boundary layer. Since the 
method of the present work is conceptually analogous to a central-difference type approach, the artificial 
dissipation must be added explicitly to the solution procedure. This is accomplished by adding harmonic 
and biharmonic terms to the governing equations, corresponding to second and fourth differences of the 
conserved variables, respectively. These dissipation terms are defined by 

D = V (t (2) a) VQ - V 2 (c (4 ) a) V 2 Q (5) 

where A is the local maximum eigenvalue of the governing equations, and and are local dissipation 
coefficients that are formulated similar to those of Jameson. 1 Furthermore, the above treatment of the 
artificial dissipation constitutes an isotropic dissipation model (independent of coordinate direction) which 
generally is only applicable to the Euler equations. For the Navicr-Stokes equations, an anisotropic model 
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is required due in part to the close spacing of points normal to the surface relative to the tangential 
distribution of points (analogous to high aspect ratio cells in structured or unstructured grid methods). 
Thus an anisotropic dissipation model was developed for use when solving the Navier-Stokes equations 
on clouds of points. 

Temporal Discretization 

Time Integration 

The governing flow equations are integrated numerically in time using an explicit multi-stage Runge- 
Kutta time-stepping scheme. 1 Typically a four-stage scheme is used to solve the Euler equations with 
the artificial dissipation evaluated only during the first stage. A five-stage scheme is used to solve the 
Navier-Stokes equations with the artificial dissipation evaluated during the first, third, and fifth stages. 

Residual Smoothing 

The Runge-Kutta time-integration scheme described in the previous section has a step size that is 
limited by the Courant-Friedricks-Lewy (CFL) condition corresponding to CFL numbers of approximately 
2.8 and 3.6 for the four-stage and five-stage schemes, respectively. To accelerate convergence to steady 
state, the CFL number may be increased by averaging the residual R with values at neighboring points. 1 
This is accomplished by replacing R by the smoothed residual R given by 

R - (V 2 R = R (6) 

where c is a constant which controls the amount of smoothing and V 2 is a harmonic operator similar to that 
used in the dissipation model. Also similar to the dissipation model, an anisotropic form of the harmonic 
operator is used when solving the Navier-Stokes equations. Equation (6) is solved approximately using 
several Jacobi iterations. Convergence to steady state is further accelerated using enthalpy damping (only 
for the Euler equations) and local time stepping. 

Boundary Conditions 

To impose the boundary conditions along the surface of the geometry being considered, ghost points 
that are located inside of the geometry are used. The locations of these ghost points are determined by 
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a simple reflection of the flow field points that are close to the surface about the edges that define the 
boundary. A similar procedure is used near the outer boundary to determine the locations of ghost points 
at which to impose the far-field boundary conditions. 

Along solid surfaces, the velocity components at the ghost points are determined from the values at 
the corresponding flow field point adjacent to the surface. When solving the Euler equations, the velocity 
components at the ghost points are determined by imposing a flow tangency or slip condition which 
requires that the velocity normal to the surface vanishes. When solving the Navier-Stokes equations, the 
velocity components at the ghost points are determined by imposing a no-slip condition which simply 
changes the sign of the values of the components at the adjacent flow field points. In either case (Euler 
or Navier-Stokes), pressure and density at the ghost points are set equal to the values at the adjacent 
flow field points. Additional conditions are imposed using the ghost points to accurately treat the shear 
stresses and heat flux terms, as well as the artificial dissipation terms. 

In the far field a characteristic analysis based on Riemann invariants is used to determine the values 
of the inviscid flow variables at the ghost points that are located outside of the outer boundary. This 
analysis correctly accounts for wave propagation in the far field which is important for rapid convergence 
to steady state. Values of the viscous flow quantities at these ghost points are set equal to the values at 
the corresponding flow field points adjacent to the outer boundary. 

Results and Discussion 

Calculations were performed first with the Euler equations and then with the Navier-Stokes equations, 
to assess the feasibility of the gridless CFD concept. The results were obtained for standard cases to 
determine the accuracy and efficiency of the methodology. All of the results were obtained on the 
Cray-YMP computer (Reynolds) at the Numerical Aerodynamic Simulation Facility located at the NASA 
Ames Research Center. 

Euler Results 

Results were obtained first by solving the Euler equations for flows about the NACA 0012 airfoil. 
The field of points that was used to model the flow about the airfoil is plotted in Fig. 2. For convenience, 


8 


the locations of these points were determined by using the cell centers from the unstructured gnd of 
Fig. 1(b), and the cloud of points for each point was taken to be the cell centers of the three triangles 
that share edges with a given triangle. To more clearly demonstrate this. Fig. 3(a) shows a close-up view 
of the unstructured grid near the airfoil nose, and Fig. 3(b) shows the points determined from the cell 
centers. Figure 3(b) also shows ghost points that are located inside of the airfoil in order to impose the 
surface boundary conditions. The computational domain has a total of 6,500 points, 134 of which are 
ghost points. It is emphasized that the unstructured grid of Fig. 1(b) was used to determine the field 
of points of Fig. 2 only for convenience. In general, any method to determine the points is acceptable. 
Efficient generation procedures to determine clouds of points have yet to be developed. 

Euler results were obtained using the points of Fig. 2 for four standard NACA 0012 airfoil 
cases corresponding to various combinations of freestream Mach number Moo and angle of attack a 
including: (1) Moo = 0.8, a = 0 ; (2) Moo = 0.85, a = 1 ; (3) Moo = 0-8, ot = 1.25 , and 
(4) Moo = 1.2, ct = 7°. All four cases were run using a CFL number of 5.0 with local time-stepping, 
residual smoothing, and enthalpy damping to accelerate convergence to steady state. Figure 4 shows 
the resulting convergence histories plotted as the log of the L 2 - norm of the density residual versus the 
CPU time in minutes. The convergence histories indicate that convergence to steady state is obtained 
in only several minutes of CPU time; thus, the method is reasonably efficient in comparison with ac- 
cepted runtimes of more conventional Euler methods (without multigrid acceleration). As further shown 
in Fig. 4, the slowest convergence is for case 2 ( Moo — 0.85, a = 1 ), which is because the solution 
contains two shock waves (upper and lower surfaces of the airfoil) of moderate strength. Therefore, it 
is slightly harder to converge the solution of case 2 in comparison with the solutions of the other cases. 
Figure 5 shows the corresponding pressure coefficient distributions C p versus the fractional chordlength 
x/c for the four NACA 0012 airfoil cases. The pressure distributions for cases 1, 2, and 3 indicate that 
the shock waves are sharply captured with only one interior point, which is somewhat surprising for a 
method that corresponds essentially to central differencing. The pressures for all four cases indicate that 
the generally accepted Euler solutions have been obtained, which suggests that the gridless CFD method 
is accurate as well as efficient for such applications. 
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Navier-Stokes Results 


Next, results were obtained by solving the Navier-Stokes equations first for a flat plate and then for 
the NACA 0012 airfoil. For the flat plate, a solution was obtained initially to assess the gridless Navier- 
Stokes capability by making comparisons with the exact Blasius solution. The field of points used in 
these calculations was generated from a structured mesh of grid points, that were uniformly distributed 
along the flat plate but clustered near the plate in the normal direction to resolve the boundary layer. 
The calculations were performed for - 0.5 and Re - 10,000. The resulting streamwise velocity 
component u (noimalized by the freestream value u e ), plotted versus the similarity variable ( y/x)y/Re , 
is shown in Fig. 6 at x/l = 0.233, 0.383, 0.533, 0.683, and 0.833. The gridless results, represented by 
the symbols, indicate that the similarity solution for a flat plate boundary layer is correctly obtained and 
that the solution agrees well with the Blasius solution. 

Navier-Stokes results were also obtained for a standard laminar case for the NACA 0012 airfoil 
corresponding to M«> = 0.5, a = 0*. and Re = 5000. The field of points that was used to model 
the flow about the airfoil was again determined for convenience by using the cell centers from an 
unstructured grid of triangles. A partial view of the unstructured grid is shown in Fig. 7(a) (generated 
from a structured grid of C-type topology), and the corresponding view of points for the gridless method 
is shown in Fig. 7(b). Closc-up views near the airfoil nose of the unstructured grid and the gridless 
field of points are shown in Figs. 8(a) and 8(b), respectively. The computational domain in the latter 
case has a total of 30,720 points, 608 of which arc ghost points. Navier-Stokes results were obtained 
using a CFL number of 4.0 with local time-stepping and residual smoothing to accelerate conveigence 
to steady state. Figure 9(a) shows the resulting convergence history plotted as the log of the X 2 ^wrm of 
the density residual versus the CPU time in minutes. The convergence history indicates that acceptable 
conveigence is obtained in less than one hour of CPU time which is reasonable considering that the 
method does not currently use multigrid to accelerate convergence to steady state. Figure 9(b) shows the 
corresponding pressure distribution which indicates that the generally accepted Navier-Stokes solution 
involving separated flow near the trailing edge has been obtained by using the gridless CFD method. To 
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more clearly see the flow solution in the trailing-edgc region, velocity vectors are presented in Fig. 10. 
The flow separates near 82% chord along the upper and lower surfaces of the airfoil, and the velocity 
vectors indicate that there are small recirculation bubbles downstream of the trailing edge. This solution 
is consistent with the Navier-Stokes solutions reported by other researchers obtained for this case using 
structured (Ref. 10) and unstructured (Ref. 11) grids. 

Concluding Remarks 

The development of a gridless CFD method for the solution of the two-dimensional Euler and Navier- 
Stokes equations was described. The method uses only clouds of points and does not require that the 
points be connected to form a grid as is necessary in conventional CFD algorithms. The gridless CFD 
approach appears to resolve the problems and inefficiencies encountered with structured or unstructured 
grid methods and, consequently, offers the greatest potential for accurately and efficiently solving viscous 
flows about complex aircraft configurations. The method was described in detail and calculations for 
standard cases were presented to assess the accuracy and efficiency of the capability. The capability 
was tested for the solution of the Euler equations and for the solution of the laminar Navier-Stokes 
equations. These solutions were found to be accurate and efficient in comparison with solutions from 
conventional CFD methods. 

Future Work 

The three-dimensional version of the gridless algorithm has been developed for the solution of the 
Euler and Navier-Stokes equations and is currently being evaluated for three-dimensional applications. An 
Euler case that is being considered is a transonic flow about the Boeing 747 transport configuration. For 
convenience, a field of points has been created from an existing unstructured mesh of tetrahedra for the 
747, using the cell centers to locate the points for use by the gridless method. The computational domain 
contains 109,805 points, 8,330 of which are ghost points. The ghost points that arc being used to model 
the surface of the 747 are shown in Fig. 11. These ghost points clearly show that the geometry includes 
the fuselage, the wings, horizontal and vertical tails, underwing pylons, and flow-through engine nacelles. 
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Fig. 1 Partial view of meshes about the NACA 0012 airfoil. 
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Fig. 2 Partial view of field of points about the NACA 0012 airfoil. 
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Fig. 6 Streamwise velocity distribution from the gridless Navier-Stokes 
solution for a flat plate at = 0.5 and Re = 10,000. 




20 


Fig. 7 Partial view of computational domains for the NACA 0012 airfoil. 
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(a) unstructured mesh of triangles. (b) corresponding field of points. 

Fig. 8 Close-up view near the nose of the NACA 0012 airfoil. 
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(a) convergence history. (b) pressure distribution. 

Fig. 9 Gridless Navier-Stokes solution for the NACA 0012 airfoil at Moo — 0-5, o = 0 , and Re — 5000. 
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Fig. 11 Ghost points for the Boeing 747 transport configuration. 
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